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The paper is devoted to the magnetic properties of isolated mesoscopic grains. We demonstrate 
that under very general conditions the electron - electron interactions in such grains can be taken 
into account by a simple interaction Hamiltonian. This Hamiltonian involves only three coupling 
constants which correspond to charging, exchange interaction, and superconducting correlations. 
The most important condition for such a description is that Thouless conductance of each grain is 
large. Under this condition sample-to-sample fluctuations of the coupling constants can be neglected. 
However, the thermodynamic properties can still remain sample-specific due to the one-electron part 
of the Hamiltonian. If the grain is made from a material which is close to the threshold of ferromag- 
netic instability the mesoscopic fluctuations of the magnetization are especially strong. Moreover, 
the situation becomes multi-stable: free energy of each grain as a function of the magnetization is 
characterized by a large number of local minima. We analyze the statistics of these minima and 
show that it possesses simple scaling properties. Numerical simulations confirm this scaling. 

PACS numbers: 73.23.-b, 71.24.+q, 75.20.Hr 



I. INTRODUCTION 

The simplest model of ferromagnetism in metallic systems was proposed by StonerQ. The magnetic ordering in this 
model takes place when the increase in the orbital energy due to the promotion of electrons to higher energy states is 
smaller than the energy gain due to the exchange interaction. As soon as this happens, the system becomes unstable 
with respect to the transition to a state with a non-zero total spin 5 and, hence, broken T-invariance. In a bulk 
system only macroscopically large 5 matters. Therefore, the onset of the instability is well determined - both the 
exchange and the orbital energies are self-averaging quantities. 

This might not be the case for a small (mesoscopic) metallic grain since both orbital and exchange energies are 
sample specific. Therefore, one should expect strong mesoscopic fluctuations of the magnetization of the grain made 
from a material that is close to the Stoner instability. This paper is an attempt to develop a theoretical description 
of these mesoscopic fluctuations. 

Consider a mesoscopic grain with an even number of electrons at zero temperature, T=0. Restricting ourselves to 
the Hartree-Fock approximation, we can speak about orbital states with energies £ a . For a weak exchange interaction 
all of the orbitals below the Fermi level Ep are double-occupied (2-orbitals), while those with £ a < Ep are empty 
(0-orbitals). This state is a singlet, 5 = 0. Let a=0,l label correspondingly the highest 2-orbital and the lowest 
0-orbital in the singlet ground state, see Fig. |l| 

For the stronger exchange a state with Stot ^ can become a ground state. Indeed, let us compare the lowest 
energies E (S) of the singlet (5 = 0) and triplet (5=1) states (see Fig. |l|). The lowest triplet state has two single 
occupied orbitals: a = 0, 1 (1-orbitals), while a > 1 and a < label 0-orbitals and 2-orbitals respectively. For the 
ferromagnetic sign of the spin exchange the two electrons on the 1-orbitals have parallel spins. Assuming SU(2) 
symmetry in the spin space of the system we can write the exchange energy of two electrons, which occupy the orbital 
states a and j3, in the usual Heisenberg form 

= -2J a phh, (1.1) 
where Si,2 are the spin operators. From this point on the superscript H indicates the quantities belonging to the 



systems with the Heisenberg form of the exchange interaction, Eq. (1.1). We can write the following expression for 



the energy difference between the lowest energy states with 5 = and 5=1, 



E?(l) - Eg(0) =£ 1 -£ - i J 01 - | J 00 . (1.2) 



1 



[In derivation of Eq. (1.2) we used the fact that 2siS2 = 5(5 + 1) — 3/2, when we add spins of two S — 1/2 particles]. 
Therefore, the triplet state becomes energetically more favorable than the singlet one, E (l) < E o (0), provided that 
J oi + 3 Joo > 2(£i — So). 

This is a sufficient rather than a necessary condition for 5 to be nonzero in the ground state. Even for Eq(1) > Eq(0), 
the spin S of the ground state is not necessarily zero. For the lowest energy among the states with a total spin 5 we 
have 



a=l 



E Q H {S)-E"(0)-'^[£ a -£i- a --J 1 - a ,i- a \- ^ T^J' 



'a/3- 



(1.3) 



1-S<a</3<S 



Here a < 1 — 5 labels 2-orbitals, 1 — 5 < a < 5 labels 1-orbitals and 0-orbitals have labels a > S. 

One can also consider situation with strong spin anisotropy. Let it be the easy axes (Ising) anisotropy, so that the 
state of an electron is characterized by the z-component s z of its spin (s z = ±1/2), and the exchange energy equals 
to 



£ Q/3 



-2J Q ps 



/3ft z (l)*z(2)- 



(1.4) 



From this point on the superscript / indicates the quantities belonging to the systems with the Ising form of the 
exchange interaction (Eq. (1.4)). The energy difference between states with zero and finite total spin in this Ising 
case is 



i i S f 1 



E 



l-S<a,P<S 
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(1.5) 



In a g ener al case the orbital energies £ a and the exchange energies J Q( g are random. In the limit 5 — ► oo summation 
in Eqs. (1.3), (1.5) leads to the self-averaging. Calling 8% = (£ a +l — £a) mean level spacing and J = (J a g) mean 
exchange coupling constant (here (...) stands for the ensemble averaging), we obtain from Eqs. (1.3) and (1.5) 



E»(S)-E?(0) = S 2 5 1 - 



3S* 



25(25-1) 



J = S 2 S 1 - JS{S +1); 



E*(S) - ^(O) = S% 



s- 



25(25- 1) 



J = S 2 6! - JS 2 



(1.6a) 
(1.6b) 



where 5 is an integer. 

Note that for a grain with an odd number of electrons the Eqs. (1.6a) and ( 1.6b| ) are slightly modified. The lowest 
possible spin is now 5 = 1/2 and the energy of a state with a spin 5 is given by 



E»(S) E» (I) = S 2 5 1 JS(S + 1) - U + ~ J; 



(1.7a) 
(1.7b) 



where 5 is a half integer. One can see that for both Heisenberg and Ising exchange interactions, the system is unstable 
at 



E'oiS) El (I) = 5^ - JS 2 - ifc + ]j, 



Si < J. 



(1.8) 



Indeed, under this condition E(S) in Eq. ( 1.6b| ) tends to —oo as 5 — > oo. Equation (1.8) is nothing but the familiar 
Stoner criterium of this instability. As soon as the parameters of the system surpass the instability threshold, the 
ground state acquires magnetization proportional to the volume of the system. 

The situation in finite systems is somewhat different. We start our discussion with a simple but instructive example 
of a grain with an even number of electrons without any disorder, i. e. £ a +\ — £ a — 8\ and J a p = J. It turns out 
that for the Heisenberg case the spin of the ground state S g is finite already at 



Si < 2J, 



(1.9) 



as can be directly seen from Eq. (1.2). Therefore, 5 ca n be greater than zero even for the exchange which is 
below the critical value given by the Stoner criterion ( |l.8D (<5i/2 = J < J < J c = Si). In this parameter domain 
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Sg = J I (2<5i — 2 J). Because Si and J are inversely prop ortional to the volume of the system, the spin of the ground 
state S g does not scale with the volume. Equation ( 1.7a ) implies that for an odd number of electrons the system spin 
S g > 1/2 (though finite) in the interval 2<5i/3 = Jo < J < J c = 8\. On the contrary, in the case of Ising interaction 
S g = (1/2) for an even (odd) number of electrons as long as Si > J. 

In more realistic models, which take the randomness in both the level spacing and the exchange interaction into 
account, the magnetization at < J < Si is essentially random. Situation far from Stoner instability was discussed 
in Refs. [2j|| In this paper we propose a theory of th e mag netiz ation near this critical point, 5i — J -C 6%. 

For finite S the functions Eo(S) from Eqs. (1.6a) and (1.6b) are random. Their statistics are determined by the 
probability distributions of E a and J a p. Below we derive the correlation function (E (Si)E (S2)) for large but finite 
Si. 2 in a realistic model of a weakly disordered or chaotic metallic grain. This correlation function is sufficient to 
describe the statistical mechanics of the grain. At high temperatures the randomness can be treated pcrturbatively 
in powers of 1/T (high temperature expansion). The low-temperature behavior of those systems is determined by the 
position and the structure of the deepest minimum of the _E(S*)-function. Therefore, one has to develop a description 
of the spin S g of the ground state: E {S g ) = mm{E (S)}. Since S g is random, one might be interested in the statistics 
of S g and of E(S) for S close to S g . The problem of evaluation of statistics of S g is not trivial. In this paper we do 
not solve it completely. We restrict ourselves to scaling analysis and numericaLsimulation. Our attempt to .construct 
an analytical description based on the Replicas Symmetry Breaking paradigmu will be published elsewhereQ. 

The remainder of this paper is organized in the following way. In section || we discuss a model of interacting electrons 
in a weakly disordered metallic grain and the statistics of E a and J a p, that follow from this model. Section III is 



devoted to the derivation of the correlation function (/^('S'O-Eo^))- The scaling analysis of the structure of the 
minima of E (S) is performed in Sec. IV. In Sec. M the results are compared with numerical simulations. Our 
findings are summarized in the Conclusion. 



II. ELECTRON-ELECTRON INTERACTIONS IN ISOLATED METALLIC GRAINS 



Let us discuss what the energy scales that determine properties of a finite system of electrons are. The single 
electron spectrum {£ a } is characterized by the mean level spacing 

Si = {£ a+ i-£ a ), (2.1) 

where (...) stands for the ensemble averaging. Another relevant energy scale of the problem is the Thouless energy 
Et ~ fr/terg, where t er g is the time it takes for a classical counterpart of an electron to cover the energy shell in the 
single-particle phase space. For diffusive and ballistic systems t erg equals to L 2 /D and L/vf respectively. Here L is 
the size of the system, vf denotes the electronic Fermi velocity, and D is the diffusion coefficient. 
The important characteristic of the system is the ratio of these two energy scales 



g = E T /8i 



(2.2) 



which is known as the dimcnsionlcss conductance. Here we discuss only metallic grains where all single-electron states 
are extended, and, thus, the dimensionless conductance is large. 



.9> I- 



(2.3) 



It is well knowriTB that in this regime the statistics of the spectrum {£ a } on the scales smaller than the Thouless 
energy Et are well described by Random Matrix Theory (RMT)EiiHl. RMT gives a quantitative description of the 
phenomenon of the level repulsion. For an ensemble of N x N, N —> oo matrices with_random and independent matrix 
elements the probability density of a realization of the spectrum {8^} is given byoO: 



P{{£y}) oc exp 



2 ^ 



\£»-£A 
Si 



(2.4) 



The parameter (3 in Eq. ( |2.4| ) can be equal to 1, 2 or 4 for orthogonal, unitary and simplectic ensembles respectively. 
The orthogonal (unitary) RM ensemble corresponds to weakly disordered grains with preserved (violated) T-invariance 
and negligible interaction between orbital and spin degrees of freedom. In what follows, we restrict ourselves by these 
two ensembles, i.e., do not consider grains with spin-orbit interaction. The latter would correspond to the simplectic 
ensemble, (3 = 4. 
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Now let us take the effects of electron-electron interaction in the grain into consideration. It turns out that a 
large class of such effects in a given disordered closed metallic grain can, under very general conditions, be described 
by a remarkably simple Hamiltonian with only three system dependent coupling constants. Let the dimensionless 
conductance g of the grain tend to infinity. We start with the simplest case when the electrons interact via a short- 
range potential 



H int (r) = \5xV5ir). 



(2.5) 



Here V cx L is the volume of the grain and A is a dimensionless coupling constant. The matrix element of this 
interaction in the basis of eigenstates (p a (r) of the noninteracting Hamiltonian is given by 

(2.6) 



M% = \5\V J dr^ a {r) V ;{r) V ^r)^{r). 

We first consider the situation without magnetic field (T-reversal invariance is preserved). Since no spatial sym- 
metries are assumed, the one particle orbitals are not degenerate an d the eigenfunctions of one-particle Hamiltonian 
ip^if) can be chosen to be real. The off-diagonal elements in Eq. ( |2.6[ ) are small as Si/ Ej- = 1/g because the integrand 
quickly oscillates. Thus, they can be neglected. O n th e contrary, diagonal matrix elements (a, 7, /1, v are equal 
pairwise) are much larger since the integrand in Eq. (2.6) is positive definite. Substituting the integrand by its mean 
value, e. g., 



we find 



07 



77 



(2.7) 



(2. 



Corrections to Eq. (2.8) ar e ne gligible for exactly the same reason the off-diagonal matrix elements are. As a result, 
in the limit Et — > 00, Eq. ( [2.8|) becomes exact! 

For the following discussion it is convenient to introduce operators of the number of electrons n a on the orbital a 
and the spin S a on this orbital 



(2.9) 



Here (a a<T ) creates (annihilates) an electron with a spin a on the orbital a and <rC are the Pauli matrices. 
Neglecting the off-diagonal terms in the interaction Hamiltonian reduces it to the form 



H int = ffM + HW + ff(3) = 

V „ fM^n* n n 4- M al 'al n a + M aa n^ n 

Z-tapcrtT! |_ c«7 acr a<T 70 , i a 7°'l < Ivl - r a a a/7 u 'ouT 1 "'j C r 1 'J"y<T i Ju 77 Uj aa u aa 1 Uj 



70-^70" i J 



To begin with, let us consider the first two terms in Eq. (2.10). Using Eq. (2.£) and the identity 



one can present + in a usual Hartree-Fock form: 



AC" -Im:\, ) „„„. -iM^SaS-, 



(2.10) 



(2.11) 



Now we can use the remarkable independence of the matrix elements of their indices, Eq. (2.8), and present Eq. (2.11) 
as 



1 "-7 
-h 2 - 2S 
2 



(2.12) 



Here S and n are the operators of the total spin and the total number of electrons correspondingly, 

S = ^Sa, h = y^n a . 



(2.13) 
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The third term in Eq. (2.10) can be treated in a similar fashion. One can write it as 



ff( 3 ) = ASif+T; f = ^a a] a al . (2.14) 

a 

It follows from Eqs. ( 2.12| ) and ( 2.1 4| ) that the Hamiltonian Hint can be presented in terms of three operators h, 



S, and T from Eqs. (2.13), (2.14), rather that in terms of all the operators h a , S a a nd a a ia a i with different a. In 
fact, this feature is not specific to the particular short range interaction Hamiltonian (|2.5| ), and is determined by the 
chaoticity of the eigenfunctions (p a (f). It turns out that <p a {r) are Gaussian random variables which are not correlated 
with each other. This is correct for eigenvectors of N x N, N — > oo random matricesllil: 

1 2 — (3 

(ip*(m)(p v (n)) = —S^Smn, {<Pn{m)(p v {n)) = — 6^5 mn , (2.15a) 

and (ip^(m)) = 0. It is also correct for the Berry ansatzEl for the wavefunctions in a chaotic grain 

1 2 — 3 

(M^M^)) = —S^h(rl - r 2 )\ (^(rl)^(r^)) = — —5^h(f{ - f 2 ). (2.15b) 

Here h(r) is the familiar Friedel function, 

Jd , (x) o-tt \r\ 

h(r)=T(d/2)^—; x=^, (2.16) 

where d is the spatial dimension of the grain, J is the Bessel function and is the Fermi wavelength. The char- 
acteristic scale of the decay of this function is of the order of the electron wavel ength, Xp , and for the purposes of 
convolution with any smooth function can be substituted by XpS(r). Equation (2.15b) becomes exact in the limit 



It is easy to see that the correlators of the eigenvectors (eigenfunctions) (2.15) are invariant with respect to an 
arbitrary orthogonal transformation performed over them: 

Mm) - E OwM"0; MV - E/ df i°^' Mn); (2.i7a) 

E °m m >°i'n = S mnS^ J dr^ O""' (fitf (r, r 2 ) = S,J (f a - f 2 ) . (2.17b) 

fi' ,m' p,' 

Therefore, the ensemble averaged part of the interaction Hamiltonian must also be invariant with respect to such 
transformations. [As to mesoscopic fluctuations of Hi nt , they disappear in the universal (RM) limit, g — > oo (N — * oo), 



and can be neglected under the condition (2.3).] There are only three operators, quadratic in the fermionic fields, 
which possess this invariance: 



,(Ti,cr 2 a ao' 1 f7 o'l,a'2 a a,cr2! 



a 2 2-i a, t 

T = J2 a ^d al , (2.18) 



and the quartic operators may be constructed only as second powers or products of these operators. Moreover, the 
Hamiltonian Hi nt should also be invariant with respect to SU(2) rotations in the spin space. Therefore, the spin may 

enter into the interaction Hamiltonian only through the combination S . Finally, Hi nt must be invariant with respect 
to U(l) gauge transformations. It means that H int can include the operator T only as a product T^T. We conclude 
that in the general case the limit of the infinite conductance g corresponds to 

H int = E c h 2 - J(S) 2 + XbcsT^T. (2.19) 

(Linear in n term is allowed, but it can be included into the one-particle part of the Hamiltonian.) We introduced 
three coupling constants E c , J and XbcS: that correspond to the three operators permitted by the symmetries. In 
the simple model with the short range interaction and preserved T - invariance, 8=1, the above coupling constants 
have the following form: 
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Ec= 2^1; 



J = 2X6i; Xbcs = XSi. 



(2.20) 



If the T - invariance is broken, (3 = 2, transformations (2.17) become unitary. The operator T from Eq. ( [2.18 ) is 
incompatible with such symmetry and, thus, Xbcs = 0, fo r (3 = 2. 



We would like to emphasize once again that Eq. ( 2.19 ) is the most general form of the Hamiltonian in the limit 
g — > 00. For instance, one can check that it correctly takes into account interactions of the two electrons on the same 
orbital. Indeed, it follows from Eq. (2.6) and Gaussian distribution of y a (r) [with correlation functions determined 
by Eqs. (2.15)], that the "double diagonal" matrix element Af™ is different from diagonal ones, given by Eq. (pTq): 



AC 



(4-{3)M^ = (4-f3)X6i. 



(2.21) 



to MZ 



On tig first sight it appears that the Hamiltonian ( [2.19 ) has to be supplemented by an additional term proportional 

H? oc XS ^2 a L,T a «.T a i,i"a,|- (2-22) 

a 

However, one notices that for any two distinct labels 0^7 there are three different off-diagonal elements M"^, MJ" 
and M 2Z, whereas double-diagonal ter m M " ° should be taken into account only once. As a result, the relationship 
(2.21) is exactly what is needed for Eq. ( |2.10| ) to be valid, provided that the sum in the right-hand side involves terms 
with a = 7 as well as those with a ^ 7. Therefore, the interaction between electrons on the same orbital does not 
violate the invariance under the rotat ion of the basis, which becomes exact in g — ► 00 limit. On the other hand, an 
additional term of the form Eq. ( 2. 2t\ ) explicitly violates this symmetry. For this reason we think that taking such a 
term into account, as it was done in Rcf.tj is incorrect. 

Let us briefly discuss the physical meaning of the three terms in the Hamiltonian Eq. (2.19). The last one exists 
only in the orthogonal case (f3 = 1) and leads to the superconducting instability provided that Xbcs < 0, i.e., there 
is an attraction in the Cooper channel. Superconducting correlations are suppressed by the magnetic field, and thus 
do not exist at (3 — 2. Here we are not interested in the effects of superconductivity and assume that the grain is a 
normal metal at T = 0. It means that Xbcs > 0. It is well known that the very same renormalization which leads to 
superconductivity at Xbcs > renormalizes the negative effective coupling constant in the Cooper channel to zero 
(see, e.g., Ref. |l3|). This fact enables us to simply neglect this interaction. 

The first two terms in Eq. ( 2.19| ) represent the dependence of the energy of the grain on the total number of the 
electrons inside and on the total spin respectively. They commute with each other and with the single-particle part 
of the Hamiltonian ^ Sa^a a Q>a.,o provided the grain is isolated. Therefore, all states of the grain can be classified by 
n and S. 

As long as we are in terested in spin structure of the system with a fixed number of electrons, we can disregard the 
first term in Eq. ( 2.19 ). As a result we arrive to a simple Hamiltonian 



H 



J(S) 2 



(2.23) 



Note that the only random component of the problem is the single-particle spectrum {£ a }, while the exchange J 
does not fluctuate! 

It should be noted though, that RM theory is just an approximation. It becomes exact only in the limit g — > 00. 
For finite g there are corrections to Eq. (2.19), which sometimes are of importance. However, these corrections at 
j > 1 do not bring essentially new physics to the problem of the small grain magnetization. For this reason we 
restrict ourselves to the Hamiltonian Eq. (2.19). 



Choosing the direction of the total spin of the system to coincide with the z-axis, we can express the energy of the 
system through the occupation numbers n a a : 



E 



H 



5> 



a ' L a.,a 



JS(S + 1), 



(2.24) 



where 



b ~ 2 



= \S Z 



(2.25) 



and Y^, n a.a = A. 



G 



Throughout this section we assumed that the metallic grain possesses rotational symmetry in the spin space. In 
addition, it is worthwhile to discuss a case in which this symmetry is broken and the inter action as only along one 
easy axes. To make connection with the energy Eq(S) that has been defined in Eq. (1.6b) we adopt the following 
Hamiltonian: 



- J{S z f 



(2.26) 



Using the notation of Eq. ( 2.25 ) we express the energy for the Ising case as 

E — ^ ^ ^ a j (j J S . 



(2.27) 



The energy Eq(S) which was introduced in the previous section can be obtained by minimizing E at fixed S with 
respect to the occupation numbers n a , a . In spite of the simple form of Eqs. (2.24), (2.27) and ( 2.25| ) the problem 
remains non-trivial, since the spectrum {£ a } is random. In the following sections we consider the effect of this 
randomness on the properties of the ground state. 



III. THE EFFECTIVE HAMILTONIAN 



We are interested in the spontaneous magnetization of the system at low temperatures. Therefore, in the following 
discussion we only compare the energies of the lowest lying states with a given spin. Hence, the state labeled by the 
total spin S is the state for which the first N/2 + S single electron states are occupied by spin up and the first N/2 — S 
by spin down electrons. We ignore the excited states with the same total spin in the following discussion. (It can be 
shown that taking these sta tes i nto account introduces a nonessential addition to the free energy). This enables us 
[based on the Hamiltonian ( [2.24 )] to write the following energy functional: 



E(S, = f £> ep(e)de - f ep(e)de - JS(S + £ - 1)), (3.1) 

JO J£ < 

where p(e) — J2 a $ ( e ~ ^a) is the one electron density of states and J is the strength of the exchange interaction, 
and parameter £ = 1,2 characterizes two limiting cases of spin anisotropy: £ = 2 corresponds to the isotropic spin 
orientation (Heisenberg model) with the degeneracy of the state 2S + 1, while £ = 1 describes the easy axes (Ising) 
model where states are only double degenerate. We have 2-orbitals at £ a < £ < , 1-orbitals at £ < < £ a < £ > and 
empty 0-orbitals at £ a > £ > . From the conservation of the total number of particles we obtain the following equations 
for £ > and £ < 

S = / p(e)de= / p(e)de. (3.2) 



The minimum of the energy functional (3.1) is the same as the ground state of the original system. The density of 
states can be represented as 

p(e) = l/5 1+ Sp(e). (3.3) 

We are interested in the value of the spin S ^> 1, which translates into the scale of orbital energies that is much larger 
than . The fluctuations of the density of states averaged over such an energy scale are much smaller than its mean 



value 1/Si, and it is sufficient to keep only the terms linear in 5p. Thus, we obtain from Eq. (3.2): 

fSSx 



£> = S — / deSp(e), 
Jo 



£ < = -S+ f de6p(e). (3.4) 

In what follows we adopt the notation: 

H = 25x-2J^8x. (3.5) 



7 



According to Eq. (1.8), the bulk Stoner instability emerges at p, = 0, so that the parameter p characterizes how close 
to the criticality thre sho ld the syst em is. 

Substituting Eqs. fl3.4|) into Eq. (3.1), we obtain: 



E(S,Z) = ^S 2 -(£-l)(Si-p/2)S-S 



Sp(e)de 



-S5i 



S(5i 



(6p(e) + Sp(-e))ede, 



Integrating Eq. ( |3.q ) by parts and using Eq. ( |3.5| ) to neglect p, as compared to Si we obtain 

E(S, = ^S 2 - (£ - 1W - [ SSl de f de x [Sp(ei) + Sp(-ei)] . 
1 Jo Jo 

Thus, we reduced the original problem to finding the minima of the random function 

E(S,0 = %S 2 -(Z-1)6 1 S + V(S), 



(3.6) 



(3.7) 



(3.8) 



where p is fixed (and small compared to Si), and V(S) is a random potential which is determined by the fluctuations 
of the density of states Sp(e) 



SiS 



V(S) = de dei [Sp(ei) + Sp(-ei)] 



(3.9) 



Below we refer to this problem as a Random Potential Problem (RPP). Such a problem is well defined provided that 
the corre lation function for random potential V(S) is given. One can evaluate the statistics of the potential V(S), 
Eq. Q, using RMT. 

The fluctuations of the density of states, Sp, averaged over energy intervals larger than Si in the ensembles of RM 
are Gaussian random variables. Therefore, using Eq. (3.9) we conclude that V(S) is a Gaussian random variable as 
well. The correlation function (V(S)V(S')) can be expressed through the correlator of the density of states. The 
correlation function of the Fourier components of Sp(e) is given byt£l: 



(SpkSp P ) = ^-S{ P + k), 



(3.10) 



where (3 is 1 or 2 depending on whether the Ha milt onian belo ngs t o an orthogonal or a unitary ensemble. Averaging 
the product V (Si)V (S2) with the help of Eqs. (|3.9|) and Eq. (3.1C), we obtain the correlation function: 



(Si + S 3 ) 



(V(Si)V(S 2 )) = 4±^ [-S?lnS? - S 2 \nS : 
ln(5i + 5 2 ) 2 + 



(3.H) 



(Si - S 2 ) 2 



HSi - s 2 ) 



Equations (3.8) and ( |3.11 ) define the random potential problem. In the next section we will show that the statistical 
description of this function possesses simple scaling properties. 



IV. SCALING ANALYSIS 



It follows from Eq. (B.ll), that the correlation function of the random potential is a homogeneous function of degree 
two: 

(V( 1 S 1 )V( 1 S 2 ))= 1 2 {V(Si)V(S 2 )). (4.1) 

Equation ( [4.l| ) means that the probability of the potential realization V(S) is the same as the probability of the 
potential 7^(7" S), and by no means expresses the scaling for the potential of a given realization. 

Our goal now is to demonstrate, that this property of the Gaussian random potential dictates a particular scaling 
form for all the moments of the free energy F of the system: 
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F(T,^;{V{S)}) = -Tin (s + ^ exp 



= -7Tn<^ 2 



dSS^ 



exp 




(4.2) 



where energy E(S) is given by Eq. (p.q). Note that in Eq. (4.2), which is valid in the continuous limit, we neglected 
unity in comparison with 2S. This is because we have already used several times the fact that S 1, and keeping 
this unity in the pre-exponential factor would be beyond the accuracy of the calculation. 

It follows from Eqs. ( |4.2| ) and ( |3.g| ) that for any given realization of the random potential V(S), the following 
identity holds: 



F (7T,7~V; {jV^S)}) + 7 Tln2 + 7 £Tln 7 = 7 F (T, /i; {V(S)}) + Tin 2. 



(4.3) 



According to Eq. ( fl.l| ) the probability of the potential realization V (5) is the same as the probability of the potential 
7^(7 S). Therefore, the moments of the free energy possess the following scaling property 



(F ( 7 T, 7 -V)> +7ZTn2 + 7 eZTn 7 = 7 (F (T, /*)) +Tln2 
(f» ( 7 T, 7 -V)> - 7" (f* (T, m)>c, « = 2, 3, 4, . . . , 



(4.4) 



where (. . .) stands for the ensemble averaging and subscript "c" means the irreducible average (cumulant). 

We can use the fact that the only available variable with the dimensionality of energy is 5i/y//3, see Eq. ( [3.11 ), and 
conclude that the moments of the free energy should be of the following scaling form 



Here we introduced scaling variable 



-T 



hi 



In 2 



2 2 



/3/x 



si 



([F(T,f,)] n ) c =[^-) if } (0,/3), n = 2,3,. 



(4.5) 



(4.6) 



that has the meaning of a dimensionless temperature. In Eq. (4.5), are dimensionless functions which can not be 
found from the scaling arguments alone. 

The stat istic s of the magnetization can be calculated with the help of the identity which follows straightforwardly 
from Eqs. O) and (O): 



S 2 



I™. 



(4.7) 



The bar in Eq. (4.7) and below stands for the thermodynamic average within a given realization of V(S). In a 
complete analogy with the derivation of Eq. (4.5), we obtain the following scaling behavior of the magnetization: 



<5 2 )-^ + ^4 1} {6,(3): 



( ) 



For n = 1 one finds from Eq. fl4.7| ) that 



A. 



G 



(n) 



,/3), n = 2,3,.. 



- h = -i6 2 ±\6-'fP(e,{3) 



(4.8) 



(4.9) 



There is no straightforward relation between and /W functions for n > 1. It is noteworthy that for the easy 
axis (Ising) case (£ = 1) the scaling functions fi (0,f3) and G^ {0,(3) do not depend on (3. 
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Let us discuss the asymptotic behavior of functions J 7 ? : from Eq. (4.5). We begin with the high temperature 



regime, 9 — > oo. One can expand Eq. (|4.2| ) up to the second order in the potential U(S), where 

U(S) = -(t-l)S + V(S), (4.10) 

to obtain 



"-^(f)-T^-I^Ml + Vm- - (4.11) 



Here we introduced the following notation 



(4.12) 



The fourth and fifth terms in the expansion ( [4.1 1|) are random quantities. Averaging them with the help of Eq. ( 3.11 ). 
we obtain 

fP^^ + 0(1/6); tf\e)=6^+0(l), (4.13a) 



7T\ 41n2 G 



2 r V 2 



0(1). 



^2 "A/ 



The numerical coefficients C 2 and D2 have the following meaning 



C 2 = f^(V(S)\ ); 



Their numerical values are 



C 2 = 



lh = P[^) {(V(S)\ )-(S.V(S)\ ( 



1 dx (1 - x) 2 1 



In - - 1 = — [#'(1/4) - #'(3/4)] w 0.3712, 



(4.13b) 



Do 



6 ln(V2 + 1) - 2y/2 In 2) w 0.5976 



where ^'(z) = ^ l ogr(z) is the second logarithmic derivative of the T-function. Substituting Eqs. ( 4.13a ) and 
( |4.13b| ) into Eq. (pL9|), we find 



G 2 1} (9) = 2/ 2 (1) (9, ft - V0\[^ + O ^ _ 



(4.15) 



Equations (|4.13aD , ( |4.13bD and pT5| ) are valid in the high temperature regime where the disorder only weakly 
affects the temperature fluctuations of the spin. There is no parametrically justified theoretical approach to anaJyae 
the situation at low temperatures. The popular approach tcuthe RPP is the Replica Symmetry Breaking ansatzBoilS. 
The results of such a calculation will be published elsewhereB. Here we restrict ourselves to a qualitative consideration, 
which yields the answers up to numerical coefficients. 

We employ arguments similar to those of LarkinO for the collective pinning of the vortex lattice, and of Imry and 
Malia for the random spi n sy stems. Let us first discuss the Ising, £ = 1, case. At the point S g of the global minimum 
of the energy E(S), Eq. ( |3.8|) , min{E(S)} — E(S g ) the random potential V(S g ) is of the same order of magnitude as 
the quadratic term and has the opposite sign. This condition can be written as 
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2 b s 



(4.16) 



At T — the entropy term is not importan t. He nce the spin at zero temperature is equal to S g and does not depend 
on the degeneracy of the state. Using Eq. (3.11), we find that 



(4.17) 



Therefore, the function tends to a constant independent of (3 as 9 — > 0. 

The estimate of the posit ion of the minimum, S^, for the Heisenberg case, £ 
slightly. It follows from Eq. (3.8) that 



2, differs from Eq. (4.16) only 



[xS g - 1 « sJ(V*(S g )). 



(4.18) 



As a result 



^2 



(4.19) 



Here a is a numerical constant of order unity, which cannot be determined within scaling considerations only. There- 
fore, the function G^ tends to a constant as 9 — > but in the Heisenberg case this constant depends on f3. 

According to Eq. ( |4.9| ), this means that /W also has a finite limit at T = 0. Moreover, the squared spin of the 
ground state fluctuates from sample to sample and the fluctuations are of the order of its average value (Sg). As a 
result, all of the functions f^ n \ must reach finite limits for both Heisenberg and Ising cases. 

Let us discuss the physical meaning of the low-temperature expansion using 



Gf ) (9) = Gf ) (0) + 



dG 



dO 



01 



(4.20) 



as an example. Both the constant term and the derivative have a simple interpretation in terms of the characteristics 
of the absolute minima of the random potential V(S). 

Consider a realization of the random potential V(S). The energy of the system ( |3.8D has only a finite number of 
minima and there are no symmetry reasons for degeneracies. Therefore, at low enough temperatures only the vicinity 
of the ground state spin S g determines all of the thermodynamic properties of the system. Close to the minimum the 
Hamiltonian can be approximated as 



E(S) ^ E(S g ) + ^ (S ~ S g f . 
The corresponding value of the square of the magnetization is 



(4.21) 



s 2 = s 2 a 



fir 



(4.22) 



Naturally, S 2 is a sa mple de pend ent quantity as well as /i r and S g . Averaging ( 4.22 ) over the realization and comparing 
the result with Eqs. (O), fl4.20|), we find 



(s g 2 ) 
\fi r 



Pi* 

dG[ 



2^1 (0): 



d6 



(4.23) 



This means that the low-temperature expansion determines the averaged location of the absolute minima and the 
curvature in this minima. 
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V. NUMERICAL SIMULATIONS 



The RPP (3.8) with the Gaussian random potential V(S), characterized by its correlator (3.11), is easily accessible 
for numerical simulations. Indeed, instead of generating the ensemble of potentials V(S) with the given correlator one 
can use the connection of V(S) with the spectra of random matrices. [See Sec. KB]. 

We carried out numerical simulations for both the orthogonal ((3 = 1) and unitary ((3 = 2) ensembles. In the 
former case we generated 1200 x 1200 symmetric (H ai — H ya ) matrices with real matrix elements H aj . These matrix 
elements were independent Gaussian random numbers with the following moments: 



= 1. 



(5.1) 



To obtain a matrix from the unitary ensemble we generated simultaneously a symmetric KcH aj and an antisym- 



metric Im_ff a7 real matrices (ReiJ Q7 = ReiJ 7Q ; lmH a 



-hxiHja) with the same dispersion as above (Eq. fly])) 



The combination ReiJ Q7 + TmiJ Q7 is a matrix element of the Hamiltonian from the unitary ensemble. 

In our analytic calculations we assumed that the mean level spacing 5\ does not depend on the location of the 
energy strip where S\ is calculated. Strictly speaking, this is not the case for the Gaussian ensembles of random 
matrices. It is well known [ll], that the density of the Random Matrix eigenvalues is distributed according to the 
Wigner semi-circle law: this density is larger in the middle of the band (e close to zero) and vanishes at the band 
edges ±6q as y/^ — e 2 . Accordingly the mean level spacing depends on the energy: 



Ji(e) =<Ji(0) 1 



(5.2) 



It is also well known that Eq. ( p.2| ) is just an asymptotic law, which becomes exact in the limit N — > oo, where N is 
the rank of the matrices; eo y/N- At finite N there are corrections to Eq. ( [5.2] ), which become most pronounced 
close to the edgesM. 

Taking all of this into account we first discarded the lowest and the highest 100 states in the spectrum of each 
Random Matrix. After that, we unfolded the rest of the spectrum according to Eq. (|]^) and obtained for each matrix 
1000 eigenstates that obey local Wigner-Dyson statistics and have uniform density. We also scaled out S\(0) and 
ended up with the mean level spacing equal to unity. 

To evaluate the realization of V(S) which corresponds to a given Random Matrix one can simply sum up energies 
of the lowest 500 + S states (filled by electrons with spin up) and the lowest 500 — S states (spins down) . Subtracting 
S 2 from the resulting sum we obtain the Random Potential V(S) in units of d±. Some particular realization of the 
random potential for different (3 are presented on Fig. ^. 

Using the generated potential V(S), we calculated the free energy, F(/z,T), Eq. (4.2) and the thermodynamic 
magnetization 



£2 = _L ^5 2 (25+l)«- 1 exp^ E{S) 



T 



z = 5>sr + i) e - 1 e* P (-ffl) 



(5.3) 
(5.4) 



for a given realization and then evaluated the ensemble average of different moments of random F(/i,T) and S 2 . All 
of the data presented below are results of averaging over 12000 realizations of random matrices. 

Figure || demonstrates the scaling properties for the mean free energy derived in the preceding section, Eq. (4.5). 
We evaluated for different values of T, fi and f3, and plotted it as a function of the scaling variable 9, Eq. (4.C). 
One can see that for the Ising case the data for different values of /z and for both /3=1 and (3 = 2 collapses on a single 
curve in accord with Eq. (4.5) . The scaling function f[ l \o) at 6 3> 1 approaches its high-temperature asymptotic 



value of 7r' 2 ln2 ps 0.070 (Eq. ( 4.13a )) within statistical errors. 

Fo r the Hcisenberg case the high temperature expansion predicts that f 2 
Eq. fT3h| )) 



(1) 



/ 2 (1) (<9,/3) w l.25x/6j3 + 0.429/3 + 0.095. 



behaves at 9 ^> 1 as (compare with 



(5.5) 



The best fit lines obtained for the numerical data are described by 
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/ 2 (1) (6>,/3) w (1.182 ± 0.002) v^3+ (0.31 ±0.02)/3 + 0.12 ±0.03. 



(5.6) 



This is in a good agreement with Eq. (5.5) taking into account the fact that there should be corrections to Eq. (4.131) 
of order \h/9 and [i (since we neglected \x as compared to 8\ using the condition in Eq. ( |3.5|)). 
Figure || illustrates the behavior of the mesoscopic fluctuations of the free energy, see Eq. (4.4). 



/f (B) = ^ ((lF(T,n)} 2 ) - (F(T, fi))' 



For the Ising case at 9 3> 1 according to Eq. ( 4.13a| ) f[ 2 ' w 0.1405 • 9 . The numerical simulations give t he slope 



of the best fit line equal to 0.1358 ± 0.0002. Therefore, the results agree well. We can rewrite Eq. (4.13b) for the 
Heisenberg case as 



V9 



0.371V9 + 2 • 0.132V^+ 0(1/Vd), 



(5.7) 



i.e f^\9,(3)/\f9 is a linear function of \/~9. The numerics indeed demonstrates such a linear dependence which can 
be best fitted by 



/f(M) 
V9 



(0.353 ± 0.006)^0 + 2 • (0.160 ± 0.003) ^ + G(l/Ve~). 



(5.8) 



Once again, the agreement is quite reasonable since there are corrections 0(/i) to the coefficients. 

Figure |B| presents the numerical results for the spontaneous magnetization ( f>.3\ ). We plot the difference between 
(S 2 ) and its high-temperature asymptotic: 



(S 2 )- 



(5.9) 



as the function of the scaling variable 9. Once again, all the curves collapse in accordance with Eq. The function 



G^-p approaches the asymptotic value within statistical errors. The collapse of the data at 9 — > justifies the order of 
magnitude estimate that led to Eq. (4.17). For the Heisenberg case Eq. (4.15) predicts the following high 9 behavior 



G ( 2 1] (9) = l.25y/9j3 + 0.858/3 + 0.190 + O 



1 

71 



The best fit lines are described by 



G { 2 1) {9) = (1.171 ±0.006) 



(0.54 ± 0.01)/3 + (0.29 ± 0.03). 



(5.10) 



(5.11) 



Again the results agree up to O(fi) in the slope and O (l/\f9^ in the intercept. There are downward deviations for 

the smallest /i at high temperature. They are likely due to the finite size effects (the magnetization becomes too close 
to its maximal value S = 500 determined by the size of the RM). As a result the scaling is violated. 

As a matter of fact, the scaling is violated when /i is too small or too large. At large fj, the typical ground state spin 
becomes of the order of unity and the condition S> 1, used throughout this article, no longer holds. As [i decreases 
the magnetization becomes of the order of the system size, S « 500, whereas in making arguments about scaling we 
assumed no upper bound on the value of S. Therefore, in obtaining the numerical values for the scaling functions at 
9 — > we used the values of /i that correspond to the ground state magnetization S g from ~ 15 to ~ 150 hundreds. 



The linear interpolation of the numerical curves at 9 
asymptotics of the scaling functions G*- 1 ', Eq. (4.8) 



results in the following values for the low temperature 



G^ (0) = .256 ± .005; G^ x) (0,^ = 1) = 1.60 ± .01; G 2 1] (0, (3 = 2) = 2.65 ± .01. 
In the Ising case the slope of this function can be determined rather well 



(i) 



dG\ 



-0.7±0.1, 



(5.12) 



(5.13a) 
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whereas evaluation of this slope in the Heisenberg case requires much better statistics. From what we had it follows 
that 



de 



= 0.1 ±0.2: 



dG%\(3 = 2) 
dB 



= 0.2 ±0.3. 



(5.13b) 



Even though the values of the slope are smaller than the statistical errors we do know the behavior of S 2 at low 
temperature. The smallness of the slope just means that the change of magnetizati on squ a red wit h temperat ure i s 
only slightly different from the one predicted by the high T expansion. Using Eqs. ( 5.13a ), ( 5.13b ) and Eqs. ( 4.23 ), 
we conclude that 



(S g 2 ^ 



0=i 



(s 3 2 y = 

(1.60 ±0.01 A (S 9 2 )J =2 



(.256 ±.005)^, 



'/V 



1.33 ±0.01)4- • 



and 



1. J 0.3 ±0.1, £ = 1 
M \ 0.7 ±0.1, C = 2 



(5.14a) 
(5.14b) 

(5.15) 



Figure |] illustrates the behavior of the averaged and rescaled zero temperature magnetization /j,{ \S g \) H and its 
square fi 2 (S 2 ) H in the Heisenberg case. The magnetization squared (S 2 } H is well described by Eq. (4.19) with the 
numerical constant a w 0.6. The contribution from disorder to the averaged magnetization (\S g \) H is an order of 
magnitude smaller tha n th e magnetization itself. Without disorder (V(S) = 0) the magnetization is the same for all 
grains S g = 1/fi (Eq. fl3.8| )). The correction due to randomness is around 9% in the orthogonal (3=1 and 4% in the 
unitary (3 = 2 case. 



VI. CONCLUSIONS 



We considered manifestations of electron-electron interactions in the properties of isolated metallic grains with large 
Thouless conductance, g 1. It turned out that the interaction effects can be taken into account by a rather simple 
interaction Hamiltonian ( |2.19 ) . We then applied this description to study the mesoscopic spontaneous magnetization 
of the metallic grains whose bulk counterparts are only sligh tly below the point of Stoner instability. In this case the 
problem maps onto the random potential problem, Eq. ( |3.8| ), where the total spin of the system plays the role of the 
coordinate. The randomness is manifested by the potential V{S) and is entirely due to the fact that the one-electron 
spectrum in such a grain is sample specific. At the same time, the fluctuations of the exchange interaction constant 
can be neglected provided th at g ^> 1. 

The correlation function (3.11) of the random potential V(S) follows directly from the Wigner-Dyson spectral 
statistics and possesses a specific invariance (4.1) under scaling transformations. This invariance dictates a particular 
scaling of the ensemble averaged thermodynamic properties of the grains as well as of the higher moments of their 
mesoscopic fluctuations. Dependence of all these q uan titi es o n tem pera ture T and on the distance from the point 
of Stoner instability ji can be determined, see Eqs. (4.5), (4.8), and 



up to some functions of the dimcnsionlcss 
effective temperature 9 = (3/iT, where (3 = 1(2) corresponds to the orthogonal (unitary) Dyson ensemble. 

According to Eq. (4.8), in the Ising case the zero-temperature magnetization typically gets reduced by a factor of y2 
when the system is driven from (3 = 1 to (3 = 2. In the Heisenberg case the average zero-temperature magnetization 
is largely determined by the non-random part of the Hamiltonian ( |3.8| ) (without disorder S g = l//i). The fluctuations 
of the magnetization become suppressed by a factor of \[2 as the system goes from the orthogonal to the unitary 
ensemble. If the grains are large enough, the transition between these ensembles can be completed in magnetic fields, 
which produce still negligible Zecman splitting. As a result, an anomalously weak magnetic field would substantially 
reduce the spontaneous magnetization in the Ising case or suppress its fluctuations in the Heisenberg one. This is 
due to the well-known fact that the unitary spectra are more rigid than orthogonal ones. However, the difference 
between the average magnetization with and without disorder is about 9% for orthogonal ((3 = 1) and 4% for unitary 
((3 = 2) ensembles (see Fig. ^) in the Heisenberg case. Therefore, in the Heisenberg case a small magnetic field .should 
suppress the average magnetization by only 5%. This is a much smaller effect than the one predicted in Refill The 
discrepancy is due to erroneous choice of the model Hamiltonian in that reference (see the discussion in Section II 
after Eq. fl2~2l])). 
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Of course, the evaluation of the scaling functions lies beyond the simple analysis. In the high-temperature regime 
it is possible to develop a regular perturbative expansion. At low temperature an analytic technique based on Replica 
Symmetry breaking paradigm can be used. The corresponding calculation will be reported elsewhere^. 

In the present paper, we analyzed the low-temperature asymptotic behavior numerically. We have shown that 
these asymptotics are determined by a single absolute minimum of a random potential (not accessible by a regular 
perturbation theory) . Using those numer ical results we were able to determine the average position and curvature for 
such minima, see Eqs. ( 5.12| ) and (5.13b). 
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FIG. 1. Spin configuration for a) the lowest energy S = state, a = 0, — 1, —2, . . . correspond to 2-orbitals. a = 1, 2, . . . 
correspond to 0-orbitals; b) the lowest energy S = 1 state, a — — 1, —2, . . . correspond to 2-orbitals, a = 0, 1 correspond to 
1-orbitals, and a = 2, 3, . . . correspond to 0-orbitals. 




FIG. 2. Several realizations of the random potential V(S) for (3 = 2. 



1G 




FIG. 3. a) The contribution to the free energy from disorder (function f[ , Eq. (4.5)) in the Ising (£ = 1) case. The 
logarithmic scale is chosen to demonstrate the fact that the scaling function goes to a constant in both low and high temperature 
regimes. The solid line is the result of the high temperature expansion. The statistical errors depend on 8 only and are plotted 
in the center of the Figure, b) Function f% in the Heisenberg case (£ = 2). The scaling functions are different for the unitary 
(/3 = 2) and orthogonal ((3 = 1) e nsembles. The s/0 scale was chosen to illustrate the agreement with the high temperature 
expansion obtained in Eq. ( 1.13b|) . Solid (dotted) line represents the predicted high 6 behavior for the unitary (orthogonal) 
case. At low temperature the functions tend to constants larger than the ones predicted by the high 9 expansion. 
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FIG. 4. a) The averaged fluctuations of the free energy of the grain in the Ising case (£ = 1) rescaled according to Eq. (4.5). 
At 9 —> the scaling function approaches a positive constant. The solid line is the high 9 asymptotic behavior described by 
Eq. (1.13a). b) The averaged fluctuations of the free energy in the Heisenberg case rescaled as in a). We divided the scaling 
functio n fj 2 by \/Q~ and plotted the ratio as a function of y/O to demonstrate the agreement with the high temperature expansion 
(Eqs. ( 1.13b)). The solid (dotted) line represents high 9 asymptotic behavior for the unitary (orthogonal) case j3 = 2(1). 
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FIG. 5. a) The contribution from disorder to the averaged magnetization squared in Ising case (£ = 1) scaled by [3fi 2 (see 
Eq. (4.8) for the definition). One can see that the scaling f unctio n has constant limits in both low and high temperature 
regimes. The solid line is a high temperature expansion (Eq. (4. IE )).The statistical errors depend on 9 only and are plotted in the 
center of the Figure, b) The contribution from disorder and —JS term to the averaged magnetization squared in the Heisenberg 
case (£ = 2) scaled as in a). The scaling functions are different for the unitary (/3 = 2) and orthogonal (/ 3 = 1 ) ensembles. The 
yd scale was chosen to illustrate the agreement with the high temperature expansion obtained in Eq. (4.15). Solid (dotted) 
line represents result of high 9 expansions for the unitary (orthogonal) case. At low temperature the functions go to constants 
higher than the ones predicted by the high 9 expansions. 




FIG. 6. □ (O) represents {fi 2 S 2 } H for /3 = 1(2). O (*) represents (n\S g \) H for /3 = 1(2). All data is for the Heisenberg 
(£ = 2) case. Without disorder the relation {fji 2 Sg} = {fi\S g \} = 1 holds. (Si is set to 1). 
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